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PREFACE 

The  work  reported  herein  was  conducted  by  the  author  as  principal  investigator  on 
subcontracts  CFSI  82-18  and  CFSI  83-7  between  Iowa  State  University  and  Calspan 
Corporation/ AEDC  Division,  operating  contractor  for  the  Aerospace  Flight  Dynamics 
testing  effort  at  Arnold  Engineering  Development  Center,  Air  Force  Systems  Command, 
Arnold  Air  Force  Station,  Tennessee.  The  work  is  a  portion  of  the  Supersonic  and 
Hypersonic  Computational  Fluid  Dynamics  project.  Project  Number  D206VW  (V32A-B3). 
The  Air  Force  Project  Manager  was  Dr.  Keith  L.  Kushman,  AEDC/DOT.  The  manuscript 
was  submitted  for  publication  on  April  30,  1985. 

The  author  wishes  to  acknowledge  the  efforts  of  Mr.  John  O.  levalts,  a  graduate  student 
at  Iowa  State  University,  whose  efforts  on  this  work  were  invaluable. 
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1.0  INTRODUCTION 

Over  the  past  decade,  sufficient  advances  have  been  made  in  the  field  of  computational 
fluid  dynamics  (CFD)  to  provide  an  extremely  useful  tool  to  engineering  practitioners.  This 
tool  is  the  CFD  technology  which  currently  sees  extensive  use  in  a  large  variety  of  aerospace 
engineering  design  and  analysis.  The  present  report  documents  an  application  of  this 
technology  resulting  in  a  new  analysis  tool  (i.e.,  computer  code)  which  is  immediately 
applicable  to  the  study  of  flow  past  specific  components  of  modern  aerospace  vehicles.  The 
present  CFD  development  effort  supports  testing  as  well  as  analysis  and  evaluation  of  re¬ 
entry  vehicles  and  high-speed  aircraft.  The  specific  technical  problem  addressed  herein 
concerns  viscous  flow  over  a  two-dimensional  or  axisymmetric  body  with  possible 
streamwise  separation. 

1.1  MOTIVATION 

Consider  the  wall/ ramp  arrangement  shown  in  Fig.  la.  This  two-dimensional 
configuration  may  represent  a  simple  model  of  a  vehicle  control  surface.  A  ramp-induced 
shock-boundary-layer  interaction  results  when  this  model  encounters  a 
supersonic/hypersonic  airstream.  The  flow  separation  resulting  from  this  type  of  interaction 
can  have  devastating  effects  on  the  control  surface  performance.  Therefore,  it  is  of  prime 
importance  to  be  able  to  predict  flows  of  this  type  accurately.  Figure  lb  illustrates  a 
cylinder/ flare  configuration.  Configurations  of  this  type  exist  and  also  are  subject  to  shock- 
induced  separation  problems.  The  present  effort  was  motivated  by  an  earlier,  unpublished 
computational  analysis  at  AEDC  of  ramp-induced  separation  on  an  actual  re-entry  vehicle 
control  surface.  The  available  Navier-Stokes  solvers  were  found  to  underpredict  the  size  of 
the  separation  zone  by  about  a  factor  of  two,  even  in  laminar  flow  where  the  turbulence 
model  is  not  involved.  Since  the  separation  zone  is  the  primary  phenomenon  degrading  RV 
flap  performance  at  flight  conditions,  and  since  high  Reynolds  number  phenomena  are  not 
adequately  simulated  in  wind  tunnels,  it  was  deemed  necessary  to  mount  an  assault  on  the 
problem  with  state-of-the-art  CFD  technology. 

The  result  of  the  effort  reported  herein  is  a  computer  code  designed  to  compute  the 
viscous  flow  fleld  on  and  about  wall/ramp  and  cylinder/flare  configurations.  This  code  uses 
state-of-the-art  CFD  technology.  To  date  the  code  has  been  tested  against  known  solutions 
for  Couette  and  Blasius  flows  and  it  has  successfully  reproduced  the  classical  solutions. 
Preliminary  solutions,  not  presented  herein,  have  been  obtained  at  AEDC  for  compressible, 
separated  ramp  flows. 
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1.2  BACKGROUND 

Numerical  computations  to  date  have  primarily  been  limited  to  the  two-dimensional 
wall/ramp  configurations.  For  example.  Refs.  1  -  4  report  numerical  computations  on 
problems  of  this  type  for  both  laminar  and  turbulent  flows.  Since  this  numerical  work  was 
reported,  advances  have  been  realized  in  the  treatment  of  the  inviscid  terms  of  the 
descriptive  equations.  These  advances  provide  for  integration  algorithms  with  substantially 
improved  stability  properties  without  the  use  of  artificial  viscosity.  These  more  robust 
algorithms  are  capable  of  solving  a  tougher  class  of  problems  than  previously  possible.  The 
remainder  of  this  report  documents  the  development  of  a  pair  of  computer  codes  utilizing 
the  new  methods. 


b.  Hollow  cylinder-flare  configuration 
Figure  1.  Specific  conflgurations  of  interest. 
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2.0  DESCRIPTIVE  EQUATIONS 
2.1  CARTESIAN  COORDINATES 


The  two-dimensional  Navier-Stokes  equations  expressed  in  Cartesian  coordinates  are 
given  by 


9x  9y 


=  0 


(1) 


where  e  =  ej  -  ev  and  f  =  f;  -  fy.  0  is  the  dependent  variable  vector  and  e,  f  are  the 
associated  flux  vectors.  The  bar  indicates  that  these  variables  are  defined  with  the  use  of  the 
Cartesian  velocity  vector  q  =  (u,  v).  The  inviscid  and  viscous  contributions  are  indicated  by 
a  subscript  i  and  v,  respectively.  Consult  Ref.  5  for  a  precise  definition  of  the  vectors  e,  f 
in  terras  of  the  primitive  flow  variables.  Equation  (1)  is  valid  only  for  2-D  flows. 
Axisymmetric  flows  are  modeled  by 


+ 

at  ax  ay 


h  =  0 


(2) 


where  e  =  Cj  -  Cy,  f  =  f;  -  fy,  and  h  =  hy.  Now  ^  is  the  dependent  vvariable  vector  and 
f,  g  are  the  associated  flux  vectors.  The  source  term,  h,  is  inversely  proportional  to  y,  the 
radial  coordinate.  This  term  can  easily  be  switched  on  for  axisymmetric  problems  and  off 
for  2-D  problems.  Note  that  the  vectors  e,  f,  h  are  defined  in  terms  of  the  cylindrical 
velocity  vector  q  =  (u,  v).  This  results  in  the  fact  that  <f>,  e,  f  appear  identical  to  e,  f  when 
u,  V  are  replaced  by  u,  v.  It  must  also  be  recognized  that  7  •  q  appearing  in  the  definitions 
for  e,  f  (see  Ref.  5)  is  given  by 


V 


(3) 


in  contrast  to  7  •  q  for  e,  f  given  by 


7 


dv 

dy 


2.2  GENERALIZED  COORDINATES 

An  attempt  has  been  made  to  construct  as  general  a  2-D/axisymmetric  Navier-Stokes 
solver  as  possible  while  maintaining  a  substantial  degree  of  operational  simplicity.  The  use 
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of  generalized  coordinates  has  contributed  to  this  effort.  Equation  (2)  is  transformed  from 
the  (t.x,y)  domain  to  the  generalized  (r,  i;)  computational  domain  to  yield 


90  ,  ►  9^  ,  j.  3e  3f  j.  -  j.  j.  h  o 

dr  3^  *  3f  3{  ‘  dr,  dr,  ^  dr, 


where  the  following  functional  forms  have  been  assumed: 


T  =  t 


(4) 


f  =  f(t,x,y) 


ij  =  >>(t,x,y) 


This  allows  for  the  possibility  of  utilizing  adaptive  grids  if  so  desired.  The  vectors  0,  e,,  fj,  h;, 
evi  fv>  hv  are  now  given: 


^  =  (e.  eu.  ev,  ee)^ 


Cj  =  (pu,  p  +  eu^.  euv,  (p  +  oe)u)^ 
fi  =  {pv,  GUV,  p  +  QV\  (p  +  pe)v)''‘ 


hi  =  y  (pv,  puv.  gv2,  (p  +  Qe)v)'^ 


Cv  ~ 


X(V  •  q)  +  Ut  +  Vx  U,) 
f^^X  +  IfxV,  +  +  TJy  U,) 


X(V  •  q)  u  +  2;<u  (f^u^  + 

"  +  +  Vx'fr,  +  £yUf 


(5a) 

(5b) 

(5c) 

(5d) 


\ 


+  TjyU,)  +  k(^xTj  +  tjxT,)  y 


(5e) 
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f 


fv  = 


M(iyUf  +  IfyU,  +  JxVj  + 

X(V  •  q)  +  2/i(fyV{  +  ijyV,) 

X(V  •  q)  V  +  flU  (fyU^  +  IJyU,  +  +  JfxV,) 


V 


+  2nV({yVj  +  I?yV,j)  +  k{fyT^  +  ^JyT. 

(5f) 


J 


0 


hv  =  - 


M(fyU{  +  IJyU,  + 

2/i(^y  V{  +  ijy  V,)  -  Ifl  v/y 


X(V  .  q)  V  +  |iu  (5yUj  +  »iyU,  +  txVj  +  ijxV,) 

+  2MV(jyVf  +  TJyV,) 


A 

/ 

+  k($yTf  +  »?yT,)y 


(5g) 


2.3  CONFIGURATIONS  OF  SPECIFIC  INTEREST 

The  flow  problems  of  specific  interest  are  those  depicted  in  Fig.  1.  Figure  la  illustrates  a 
wall/ramp  intersection  encountering  a  high-speed  viscous  flow.  The  possible  separation 
resulting  from  this  complex  interaction  makes  this  problem  of  more  than  academic  interest 
for  aerospace  vehicles  in  which  control  surfaces  are  used.  Figure  lb  illustrates  a 
cylinder/flare  conflguration  which  exhibits  similar  separation  characteristics  but  under 
different  conditions  due  to  the  flow’s  radial  relief.  The  current  approach  is  to  identify  a 
computational  region  on  these  configurations  and  solve  Eq.  (4)  on  this  domain  subject  to 
proper  boundary  conditions. 

2.4  BOUNDARY  CONDITIONS 

In  either  the  wall/ramp  case  or  the  cylinder/flare  case  the  computational  domain  is 
illustrated  in  Fig.  2.  The  boundaries  of  this  domain  are  labeled  A,B,C,D  corresponding  to 
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Figure  2.  ComputationRl  domain  layout. 


body,  outflow,  upper,  and  inflow,  respectively.  Boundary  conditions  consistent  with  the 
viscous  flow  requirement  are  utilized  on  boundary  A.  That  is,  the  no-slip  condition  requires 
u  =  V  =  0  on  this  boundary.  The  surface  temperature  is  either  prescribed  on  this  boundary 
or  is  determined  from  an  adiabatic  wall  condition  by  the  equation 

a  ^  +  (I  -  S)(T  -  T«)  =  0 

where  T*  is  the  prescribed  wall  temperature  value  and  T  is  the  calculated  wall  temperature,  a 
Is  0  for  the  fixed  wall  case  and  1  for  the  adiabatic  case.  Since  n  is  the  direction  normal  to  the 
boundary  A,  it  lies  in  the  direction  of  Vjj.  The  actual  numerical  implementation  of  this 
condition  is  reserved  for  Section  3.5.  Finally,  a  condition  on  the  normal  pressure  gradient  is 
enforced  to  provide  closure  on  the  wall/ ramp  boundary.  The  simplest  condition  is  dp/dn 
=  0  which  is  consistent  with  boundary-layer  theory.  Results  illustrated  in  Ref.  5  show  this  to 
be  adequate  for  the  cases  considered.  However,  Ref.  6  shows  contrary  results  for  a  turbulent 
separated  case  near  the  reattachment  region.  For  the  present  work  in  the  interest  of 
simplicity,  the  simple  zero  normal  pressure  gradient  is  enforced.  For  the  more  general  case, 
the  vector  momentum  equation  is  dotted  with  the  normal  direction  to  yield  a  rather 
complex  expression  for  dp/dn  which  must  be  solved  for  the  wall  pressure.  This  expression  is 
given  by 
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^  I  +  V,)]  +  2S,J7,-^(aC77^U,) 

"t"  2iyTJy  +  Ox^y  '*'  ’*'  ^y^if)] 

+  (^  +  i7y)  ^  [x(v  •  q)]  +  vl-^  [2/1  (^xUj  +  i;,u,)] 

+  1??^  [2/*(iyV£  +  »Jyvj] 

+  2lIx1Jy  ^  [^  (ix^t  +  iyU{  +  »;xV^  +  TfyU^] 

r 

This  completes  the  ideology  behind  the  viscous  boundary  conditions  on  boundary  A.  As 
mentioned,  numerical  implementation  of  these  is  discussed  in  Section  3.5. 

Boundary  B  is  an  outflow  boundary.  The  simplest  treatment  of  boundaries  of  this  type  is 
to  enforce  either  a  zero  streamwise  first  derivative  or  a  zero  streamwise  second  derivative  of 
all  dependent  flow  variables.  This  is  stated  as 


or 


(7) 


This  is  equivalent  to  a  zeroth  order  or  first-order  extrapolation  of  <t>  to  boundary  B  from  the 
interior  of  the  domain  along  an  1;  =  constant  curve.  It  should  be  noted  that  while  this 
condition  is  simple  and  commonly  used  with  success,  it  is  not  physically  correct.  This  point  is 
addressed  further  in  Section  3.5. 
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Boundary  C  is  generally  thought  of  as  having  fixed  dependent  variables.  However,  an 
alternate  option  is  considered  in  the  present  work.  If  flow  exits  the  computational  domain 
through  boundary  C,  then  it  is  technically  correct  to  fix  only  one  flow  variable  on  this 
boundary.  The  rest  must  be  computed  from  the  descriptive  equations.  Pressure  is  chosen  as 
the  fixed  variable  for  the  present  work. 

Boundary  D  is  an  inflow  boundary  in  which  the  dependent  variables  are  fixed.  This  is 
technically  correct  from  a  signal  propagation  viewpoint  only  in  the  supersonic  portion  of 
this  boundary.  However,  practice  has  demonstrated  this  to  be  an  adequate  boundary 
condition  for  flows  of  the  type  considered  here. 

The  numerical  aspects  of  all  of  these  boundaries  are  best  left  for  a  discretized  setting  and 
are  presented  in  Section  3.S. 


3.0  DISCRETIZATION 

3.1  GRID  GENERATION  AND  METRICS 

The  numerical  solution  of  Eq.  (4)  requires  values  for  the  metrics  J,,  jj,,  Jy,  j/y.  The 
time  metrics  are  zero  for  the  present  work  since  a  moving  grid  is  not  presently  considered. 
The  remaining  metrics  are  obtained  from  the  grid  point  locations  and  the  metric  identities. 

j  j 


“  A 


_  tii  X  rf  *  j  _  _ 


Vx 


^  r{-  X  r{  .  i  ^  _y^ 


Vy  = 


ff  X  ff  ‘  j 
J 


(8a) 


where 

J  =  rj  X  r,  .  rf  =  x^y,  -  x,yj  (8b) 

since  for  this  2-D  domain  f  =  z  and  r^  =  (0,  0, 1).  Now  since  all  calculations  are  performed 
in  the  computational  domain,  these  metrics  are  easily  computed  once  the  values  of  x^,  y^,  x,, 
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are  known.  These  values  are  determined  with  finite  difference  formulae.  The  particular 
formulae  chosen  depend  upon  how  the  flux  derivatives  are  treated.  A  discussion  of  this  is 
reserved  for  Section  3.4  but  regardless  of  the  specific  formulae  chosen,  the  grid  point 
locations  must  first  be  known. 

The  approach  used  in  the  present  work  is  to  break  the  grid  generation  problem  up  into 
four  distinct  parts: 

1 .  Boundary  point  distribution 

2.  Forcing  function  calculation 

3.  Relaxation 

4.  Redistribution 

Each  of  these  is  now  discussed.  Boundary  point  distribution  simply  means  that  the  analyst 
decides  where  to  place  the  grid  points  on  all  four  computational  boundaries  beginning  at 
boundary  A  and  moving  counterclockwise  to  boundary  D  (refer  to  Fig.  2).  The  analyst  must 
create  the  (x,y)  coordinate  values  for  each  grid  point  in  this  counterclockwise  cycle.  This  is 
typically  accomplished  by  a  separate  computer  program  which  distributes  points  on  the 
boundaries  according  to  some  measure  of  expected  flow  Held  activity.  The  points  need  only 
to  be  equally  spaced  on  boundaries  B  and  D  at  this  time  since  the  redistribution  phase 
redistributes  points  along  all  £  =  constant  curves  according  to  an  input  from  the  analyst 
regarding  boundary-layer  resolution.  If  J,K  are  used  to  represent  the  number  of  points  in  the 

n  directions,  respectively,  and  j,k  are  the  indices  of  these  points,  then  part  1  of  the  grid 
generation  process  produces 


Xj.i,  yj,!  for  j  =  1,  2 . J 

xj.k  ,  yj,k  for  k  =  1,2 . K 

Xj,K  »  yj.K  for  j  =  1,2,  ....,  J 

xi.k  .  Yl.k  for  k  =  1,2 . K 

See  Fig.  2  for  reference.  Part  2  of  the  grid  generation  process  uses  the  work  of  Ref.  7  to  force 
the  interior  points  in  the  grid  to  respond  to  the  boundary  point  distribution  created  in  part  1 . 
The  general  grid  generation  approach  used  is  due  to  Refs.  8-9  in  which  two  coupled  Poisson- 
type  partial  differential  equations  are  solved  for  the  x,y  values  at  all  interior  points  given  the 
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x,y  values  on  the  boundaries.  The  grid  equations 

=  P(i,  r,) 

=  Q(«,  ri) 

are  expressed  in  computational  space  as 


arj£  -  2j3rj,  +  7r^  +  J^Prj  +  j2Qr,  =  0  (9) 

where  r  =  (x,y),  a  =  x^+ |8  =  x^x,  +  and  7  =  x|  +  y|.  Equation  (9)  represents  two 
coupled  equations  which  are  solved  by  relaxation  once  the  functions  P,  Q  are  determined. 
These  functions,  in  general,  provide  for  arbitrary  point  clustering  but  in  the  present  work  are 
used  to  force  interior  points  to  respond  to  boundary  point  placement  as  suggested  in  Ref.  7. 
Again,  using  the  index  notation,  the  P,  Q  functions  are  defined  as 


P  = 

Q  =  iKi.  u)  |V,|' 


where  tf),  ij)  are  determined  strictly  from  boundary  point  information  as  follows. 
Along  boundary  A,  $  is  defined  as 


or 


>7=1)=- 


St 


$j.i  =  -  2(Sj+]  +  Sj_i  -  2sj)/(sj  +  i  -  Sj_i)  (10a) 

in  flnite-difference  form.  Note  that  f  =  j,  17  =  k  are  assumed  so  that  =  Aij  =  1 .  The  s  in 
these  equations  is  simply  the  arc  length  along  the  boundary.  This  is  determined  by  evaluating 
the  quadrature 

, _ 

Sj  =  J  V(*{  +  yt)j,i  d? 


Likewise  on  boundary  C, 

^j,K  =  -  2(sj+i  +  Sj_i  -  2sj)/(s3+i  -  Sj_,)  (10b) 
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where 


Sj  =  1  V(Xf  +  y|)i,Kdi 


Equations  (10a)  and  (10b)  are  applied  only  to  interior  boundary  points  (i.e.,  j  =  2,  3 . . 

J  -  1).  The  interior  field  values  of  $  at  j  for  k  =  2 . K  - 1  are  simply  found  from  the 

interpolation  formula 


♦j,l.  =  ♦i.l  +  (^)(*i,K  -  *j,l) 


for  j  -  2,  3 . J  -  1  and  k  =  2,  3 . K  —  1.  On  boundary  B, 


=  -2(sk+i  +  Sfc-i  “  2sk)/(sit+i  -  Sit_i) 


where 


I 


and  on  boundary  D, 


^l,k  —  ~2(Sk+i  +  Si[_i  “  2S|[)/(Slc+|  —  S][_]) 


where 


Sk  =  I  +  ypl.k  di? 

I 


Then  for  the  interior  grid  points 

lAj.k  =  -  ^^i,k) 

for  j  =  2,  3,  J  —  1  and  k  =  2,  3 . K  -  1.  Now,  ^  are  known  at  all  interior  grid 

points.  An  initial  grid  is  then  prescribed  to  begin  the  relaxation  process.  The  P,Q  functions 
are  computed  as  indicated  from  ^  and  V{,  Vt;.  Part  3  of  the  grid  generation  process  is 
then  carried  out  to  produce  a  converged  grid.  The  grid  resulting  from  this  process  will  not  be 
sufficiently  refined  near  the  17  =  1  boundary  to  resolve  the  viscous  features  of  the  flow  field. 
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The  redistribution  process  in  part  4  is  designed  to  resolve  this  deficiency.  Following  part  3,  a 
grid  such  as  that  shown  in  Fig.  3  is  produced.  The  points  are  now  redistributed  along  each  ^ 
=  constant  curve  to  satisfy  some  user  prescribed  near  wall  point  spacing.  Consider  the  ^ 
=  constant  curve  given  by  ^  =  4  in  Fig.  3.  The  first  requirement  of  the  redistribution 
process  is  that  the  same  total  arc  length  be  reproduced.  This  arc  length  is  computed  by 


ij=K 


Sj  =  \  V(x;  +  y^)j.k  dn 


Figure  3.  Redistribution  along  {  =  constant  lines. 


In  general,  the  desire  of  the  present  work  is  to  redistribute  points  to  satisfy  the  relation  (see 
Ref.  10) 


Sj(if)  =  Gj(7f)Sj 

where 


(Ha) 


(Hb) 


with  1)  =  (j?  -  1)/K  -  1).  This  is  to  be  subject  to  the  constraint  that  Sj(2)  -  Asj  which  is 
prescribed  by  the  analyst  as  direct  control  of  the  point  spacing  near  the  wall.  Then 
application  of  this  constraint  to  Eq.  (11)  yields 
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where 

-(t^) 

The  value  of  |3j  is  then  determined  from  this  equation  using  a  Newton-Raphson  iteration. 
The  points  along  each  {  =  constant  curve  are  then  redistributed  using  these  /3j  values  in  Eq. 
(11)  and  the  grid  generation  process  is  complete.  An  example  of  a  grid  generated  by  this 
method  is  illustrated  in  Fig.  4. 


Figure  4.  Completed  grid. 


3.2  FLUX  DIFFERENCE  SPLITTING 

The  concept  of  splitting  the  convection  terms  was  introduced  at  least  as  early  as  1952  in 
Ref.  1 1 .  Since  then  the  idea  has  been  expanded  in  a  number  of  ways  by  various  investigators 
(see  Refs.  12-19).  The  flux  difference  splitting  concept  used  in  the  present  work  closely 
resembles  that  of  Ref.  19  but  the  present  work  blends  this  concept  into  the  framework  of  the 
implicit  solution  of  the  Navier-Stokes  equations  in  factored  delta  form  and  the  splitting 
concepts  are  developed  in  a  more  systematic  and  rigorous  fashion  than  in  Ref.  19. 
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Consider  Eq.  (4)  expressed  in  the  form 


dr  V^‘  d(  3f  af 


d<t>  ^  dc:  afj 

'  d7i  *  dtj  ^  dv 


dv  ) 


= 


acy  ac  acv  af„ 

+  -57-  +  »Jx  +  »?v 


ae  af 


3v  'y  ai? 


(13) 


Let  A,  B  represent  the  Jacobian  matrices 


acj  afj 

a^  a^ 

respectively  and  let  A*  =  +  ^^A  +  and  B*  =  ij,I  +  tj^A  +  ijyB  where  I  is  the  identity 

matrix.  Then  Eq.  (13)  becomes 


w  +  A*i^  +  +  hi  = 

dr  af  aTj  '  *  af 


dfy  ,  aCy 
— ^  ’Jx  — ^ 


af 


ai7 


3fy 

drj 


+  h„ 


(14) 


For  an  inviscid  problem,  the  right-hand  side  of  Eq.  (14)  is  zero  and  the  matrices  A*,B'* 
determine  the  characteristic  structure  of  the  resulting  hyperbolic  system.  It  is  for  this  reason 
that  the  splitting  concept  in  the  present  work  is  applied  only  to  the  inviscid  terms.  The 
diffusive  terms  on  the  right  side  of  Eq,  (14)  are  not  considered  to  contribute  to  the  signal 
propagation  direction  information  in  the  fluid.  Note,  however,  that  the  viscous  contribution 
does,  in  fact,  include  terms  which  can  be  expressed  as  first  derivatives  of  <t>  multiplied  by 
various  functions  implicitly  dependent  upon  0  for  the  case  of  an  axisymmetric  problem.  For 
example,  derivatives  of  e^  and  fy  which  involve  V  •  q  contribute  due  to  the  v/y  term  existing 
in  r  •  q.  A  portion  of  the  hy  term  also  contributes.  It  is  therefore  possible  to  include  these 
contributions  in  the  matrices  A*  and  B*,  thereby  altering  the  characteristic  structure  of  the 
governing  equations.  However,  as  mentioned,  this  idea  was  avoided  in  the  present  work. 

The  development  of  the  split  equation  requires  first  that  the  system’s  characteristic 
structure  be  exploited  in  detail.  This  is  accomplished  by  diagonalizing  A*  and  B*  as  follows. 
Let  represent  the  matrix  whose  rows  are  the  left  eigenvectors  of  A*.  Let  T^  represent  the 
matrix  whose  columns  are  the  right  eigenvectors  of  A*,  and  let  represent  the  diagonal 
matrix  whose  nonzero  elements  are  the  eigenvalues  of  A*.  Likewise  let  T“ ',  T^,  represent 
left  eigenvector,  right  eigenvector,  and  diagonal  eigenvalue  matrices  of  the  matrix  B*.  The 
A*,  B*  can  be  expressed  as 
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A*  -  TfA(T,-‘,  B"  =  T,A,T 


-1 

11 


Since 


and 


dtjt  90  de:  dfj 

H  ^  If  ~df 


90  90  dcj  df[ 

=  ’Jtir  +  ’fx-sr  +  ’'y 


Bit  dtf  '  di)  y  3tj 


the  quantities  30/ d£,  30/ dr;  may  be  expressed  as 


—  -  A*~'( t  f  f 


?=B*- 

aij 


f  S* 

dCj 

8f|  \ 

dij 

But  from  Eq.  (IS), 


A*  '  = 


lx-1 


B*  =  V-^T 
Equation  (14)  requires  the  product  A*  30/3f,  which  becomes 


(15) 


(16a) 

(16b) 

(17a) 

(17b) 

(18) 


by  substitution  of  Eq.  (15).  This  term  may  now  be  split  into  a  part  involving  positive 
eigenvalues  plus  a  part  involving  negative  eigenvalues.  This  is  desired  for  the  purpose  of 
determining  which  grid  points  to  use  in  finite-differencing  the  ^  derivatives  in  order  to 
maintain  consistency  with  the  preferred  propagation  directions  in  the  hyperbolic  part  of  the 
governing  equations.  Since  is  the  diagonal  matrix  of  eigenvalues,  let  represent  the 

parts  involving  positive  and  negative  eigenvalues,  respectively.  Then 

Aj  =  +  Af 
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and  Eq.  (18)  becomes 

which  by  substitution  from  Eqs.  {16a)  and  (17a)  becomes 


de; 


or 


+ 


£!l 

dk 


with  the  help  of  the  definitions 


Q*  .  Aj-  'T(- '  Q-  =  T,A,-A,-  'Tf ' 

Now  an  Identical  treatment  of  the  B*  term  in  Eq.  (14)  with  the  definitions 
R  +  =  T,A;  a  -  'T- '  R  -  =  T,A  -  A  -  'T;  ' 

yields  the  continuous  but  split  version  of  Eq.  (14). 
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The  inviscid  counterpart  of  Eq.  (19)  is  presented  in  Ref.  19.  The  present  development  of  this 
equation  is,  however,  entirely  different  from  that  of  Ref.  19.  Clearly  for  Eq.  (19)  to  be 
consistent,  the  following  requirements  hold 

Q+  +  Q-  =  1 

R+  +  R-  =  I 

Note  that  the  construction  ofQ'*’,Q",R'^,R'  requires  an  inverse  of  diagonal  eigenvalue 
matrices.  Consider  for  example, 

Q"  =  VfAf'Tf' 


and 

Aj  =  diag  +  5,u  +  fyV,  +  {*u  +  (yV, 


+  cVfx  +  +  txU  +  ^yV  -  c^/^x  + 


where  c  is  the  sonic  velocity  of  the  fluid.  The  inverse  of  Aj  does  not  exist  if  +  ^;,u  +  (yV 
=  0  or  I  {(  +  +  fyV  1  =  c  prevent  these  circumstances  from  causing 

difficulties,  the  products  Aj+Af*  and  Aj“Aj-'  are  treated  in  the  following  way.  Let  Aj 
=  diag  (Xi,  X2.  X3,  X4)  and  let  =  Af  Af '  =  diag  (Xj+Aj,  Xj+ZXj.  XjVXj,  X^/ig.  Then 


Likewise, 


0  Xi  <  0 


XjVXj  =  V=  \  1/2  Xi  =  0  for  i  =  1, 2,3.4 
1  X;  >  0 


0  Xj  >  0 


XjVXi  =  \r=  <  1/2  Xi  =  0  for  i  =  1, 2,3,4 


I  Xi  <  0 


(20a) 
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This  way,  the  singular  character  of  the  eigenvalue  matrices  at  selected  points  is  avoided.  The 
same  procedure  is  used  to  determine  R+  and  R"  corresponding  to  the  ij  direction. 

Equation  (19)  has  some  very  special  features.  First,  the  matrices  Q'^,  Q“,  R+,  R“  act  as 
complex  weight  factors  weighting  the  various  fluid  signals  (waves)  according  to  the 
appropriate  propagation  directions.  This  is  accomplished  in  a  discretized  setting  by  using 
backward  finite-difference  approximations  for  derivatives  weighted  by  positive  eigenvalues 
and  forward  finite-differences  for  derivatives  weighted  by  negative  eigenvalues.  This  is 
consistent  with  the  wave  propagation  physics  since  positive  eigenvalues  mean  characteristic 
signal  paths  toward  the  positive  coordinate  direction  with  the  opposite  true  for  negative 
eigenvalues.  In  addition  to  this  characteristic- like  structure,  the  equation  retains  its 
conservation  law  form,  thus  providing  a  needed  shock-capturing  capability  to  the  resulting 
algorithm.  See  Ref.  20  for  a  discussion  of  the  connection  between  conservation  law  form 
and  the  computation  of  weak  solutions. 

3.3  INTEGRATION  SCHEME 

The  implicit  integration  algorithm  used  in  the  present  work  is  developed  in  this  section. 
Before  beginning  the  development,  some  comments  regarding  notation  are  in  order. 
Variables  are  generally  functions  of  (t,  17).  Thus,  when  writing  a  variable  it  is  important  to 

identify  which  value  of  r  and  at  which  point  in  space  (f,  77)  the  variable  is  applicable. 
Generally  this  is  accomplished  in  discrete  work  with  an  index  notation.  In  this  work,  an  n 
superscript  is  used  for  the  t  index  and  j,  k  subscripts  are  used  for  the  jj  indices, 
respectively.  Therefore,  the  variable  at  the  point  $  =  j  Af,  ij  =  kAj;  and  t  =  nAr  would  be 
identified  by  0]*^  where  A^,  A»j,  and  At  are  the  computational  grid  point  spacings  and  the 
time  increment,  respectively.  Because  of  the  notational  complexity  of  this  and  other 
sections,  a  shorthand  notation  is  followed.  If  the  term  0  is  encountered  with  no  indices,  it  is 
assumed  to  be  If  ^j  +  i  Is  encountered,  the  t  index  is  assumed  to  be  n  and  the  17  index  is 
assumed  to  be  k.  Simply  put,  only  indices  not  equal  to  n,  j.  k  are  explicitly  written.  This 
approach  greatly  enhances  the  readability  of  this  and  other  sections.  In  addition,  the 
symbols  A,  V  imply  first-order  forward  and  backward  differences,  respectively.  No  subscript 
implies  differencing  in  the  t  coordinate  and  i,  77  subscripts  imply  general  differencing  in  the 
77  coordinates. 

With  this  notation,  consider, the  equation 

0"'^'  =  0  -f  (1  -  i3)AT0^  +  (21) 

If  ^  =  0  an  Euler  explicit  scheme  results,  /3  =  1/2  corresponds  to  trapezoidal  time 
differencing,  and  ^  =  1  gives  Euler  implicit  differencing.  Consider  0t"'''’  from  Eq.  (19). 
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^n  +  i  =  +  Q-J  +  A0  +  0(A^)) 

+  +  AA0  +  0(A2)) 

+  +  BA0  +  0(A2))J 

+  |^R+  +  R-J  +  A<>  +  0(A2)) 


(22) 


+  ^r‘ 


y  (cj  +  AA0  +  0(A2)) 


+  BA^  +  0(A2))J  +  hf  +  ' 

[ 


- 

■  t*  a 


n  + 1  aril  ■*■  1 

+  1 
y 


+  ijr*  — 


arll+  ] 

-  +  h5+» 


] 


Note  that  Q"*"  +  Q  may  be  evaluated  at  n  or  n  +  I  since  in  either  case  the  result  is  the 
identity  matrix.  Now  with  the  additional  linearizations 


hp+*  =  hi  +  HjA^  +  0(a2) 
e;+'  =  e^  +  CA^^l  +  0(a2) 
fJ+‘  =  +  DA^  +  0(A2) 

h;+‘  =  h^  +  H^A0  +  0(a2) 


dhi 


where 
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C  = 


dCy 

d0 


D  = 


3fv 

d<^ 


and 


and  also 


Hv  =  — 


S4> 


=  fi  +  ASj 

=  ’Jr  +  AtJx 

=  ^1  +  Ajj, 

=  fy  +  Afy 

=  f X  +  A^x 

’ly'*’’  =  IJy  +  Al7y 

substitution  into  Eq.  (22)  yields 


-  «-  =  [[(q*  +  q-)(s4i  +  {,|a  +  f,iB)  +  (r^  .  R-) 

(’■i*  + ’»l®)  ■"  “] 

■  +  ”■] 


+  (q^  +  Q-)(A^t^j  +  Af,  +  (r+  +  R  ) 

“  r  A^a  Cyj  +  +  A?J,  +  Aj)y  f,,  1  +  0{a2) 
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where  alE  terms  of  0(A^)  were  ignored.  Substitution  of  this  result  into  Eq.  (21)  gives 
[l  +  +  Q-)(j,|l  +  f,lA  +  f,AB)  +  (r*  +  R-) 


(’■  i  ‘  ^  i  ^  +  ’V  ^  b)  +  H,  -  (j.  A  c  +  f,  A  D  +  ±  c 


A<f>  =  —  jSAr  [<«-.  Q“)(Aft<^j  +  A$,ei^  (23) 


+  ASyfj^)  +  +  R  )(A?it^,  +  ATf^Cj,  +  AJ7yfi,)J 

+  ^Ar^Af^e^^  +  Afyfy^  +  +  Aijyf^^)  +  0(A2) 


Equation  (23)  represents  the  general  time  differencing  formula  for  interior  points  of  the 
computational  domain.  This  scheme  is  second-order  accurate  in  t  if  /3  =  1/2.  The  spatial 
accuracy  of  the  steady-state  solution  is  dependent  only  upon  the  f,  jj  differences  used  on  the 
right-hand  side  of  Eq.  (23).  It  should  be  noted  that  only  first-order  spatial  difference 
approximations  are  required  on  the  left  side  of  Eq.  (23)  to  maintain  an  accuracy  consistent 
with  the  second-order  time  differencing.  If  the  transient  accuracy  of  the  solution  is  not  of 
primary  importance,  it  is  worthwhile  to  drop  the  source  terms  and  viscous  terms  on  the  left 
side  of  Eq.  (23).  The  result  is  first-order  accurate  in  r,  and  good  stability  properties  remain 
so  long  as  the  convective  terms  are  treated  according  to  the  splitting  philosophy.  The  result 
of  dropping  these  terms  and  factoring  the  remaining  operator  yields 


Jl  +  0ArQ+  (ftVjI  +  +  gyV^fi)  +  ^ATQ-(fiAj 

X  -H  (3AtR+(ijiV,I  -I- 
+  /3ATR”(i7tA,I  +  i;*A„A  -h  iJyA^B^l  A^ 


')! 
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-  0At 


^Q  +  ^AfiVj^  +  +  A^yVjfi^  +  +  Af^AjCj  +  AfyAffj) 


+  R  +  (Aij,V,<i)  +  Aiy.V^ej  +  ArjyV/j)  +  R~(Atj,A^0  +  Atj^A^C,  +  AijyA^fjjJ 


where  tj>r  is  given  by  solving  Eq.  (19)  lo  yield 


=  (^x5^ev  +  5yfi{fv  +  +  I7y6^f^  + 

-  [<3*(«.'’(  *  +  f.rt*ei  +  £,ri*fi) 

+  Q'(«,rt-*  +  £,rj-ei  +  £,r;-f,) 

+  R+(>},r;^  +  7,^r;ei  +  i?yr,ni) 

+  R"(’iir,7./.  +  »?,r-ei  +  iTyP-f,)  +  hj 

The  factoring  introduces  an  error  of  the  same  order  as  the  error  in  the  second-order 
differencing  scheme  which  is  appropriately  ignored.  The  operators  5j,  5,  are  difference 
operators  used  with  the  diffusive  terms  and  are  now  precisely  defined.  Consider  the  typical 
term  5j  e^.  A  quick  look  at  Eq.  (5e)  reveals  that  5j  e^  will  always  take  one  of  three  forms.  Any 
of 

6j(cg^s),  orSj(sg)  (25) 
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where  c  and  $  represent  quantities  dependent  only  upon  the  elements  of  and  g  represents  a 
geometry  related  quantity.  Then  these  three  forms  are  approximated  by  the  following 
second-order  difference  formulae: 


(cg)j+i/2(sj  +  :  “  Sj)  “  (cB)j-i/2(sj  -  Sj-i) 


(26a) 


_  (cg)j+i  (sj  +  i^lc-H  ^i  +  l.k-l)  (^8)j- l(^j- l.k+ 1 


4A(A)) 


«,(sg)  =  ,2^, 

where  half-point  evaluations  are  simply  averages  of  the  neighboring  two  points.  Similar 
formulae  are  used  for  the  6,6^  type  terms.  The  operators  r|;  r*in  Eq.  (24b)  are  one-sided 
difference  formulae  and  can  be  either  first-order  or  second-order,  depending  on  the  steady- 


state  solution  accuracy  required.  For  first-order  accurate  steady-state  solutions, 

rj+  =  Vf.  r;  =  v,,  and  r,-  =  a,  (27) 

while  achieving  second-order  accurate  steady-state  solutions  requires  the  use  of 

r,+  =  |(3V,  -  {28b) 

rr -7(3-3! 
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For  example. 


rf  *  =  -j  Me. 

In  practice,  the  value  of  0  in  Eq.  (24)  is  set  to  zero  if  the  input  value  of  jS  is  not  1/2.  This 
is  consistent  since,  if  i3  =  0  or  0  =  1,  then  the  resulting  algorithm  is  only  first-order  time 
accurate  anyway,  so  the  0  terms  are  the  same  order  as  the  error.  Therefore,  they  may  just  as 
well  be  eliminated.  If  computation  of  the  0  terms  is  required,  the  forward  time  differences  in 
Eq.  (24a)  may  be  replaced  by  backward  time  differences  without  compromising  the  accuracy 
of  the  algorithm.  It  is  sometimes  more  convenient  to  generate  differences  based  on  past 
geometric  information  rather  than  future  information.  This  is  done  in  the  present  work. 
Note,  however,  that  all  0  terms  vanish  anyway  for  a  stationary  grid. 


The  difference  scheme  represented  by  Eq.  (24)  uses  fully  general  geometric  description 
and  is  capable  of  handling  adaptive  grids.  The  laws  governing  the  grid  motion  are  not 
supplied  but  could  be  incorporated.  The  scheme  can  be  either  first-  or  second-order  accurate 
in  space  and  is  currently  first  order  in  time.  The  addition  of  the  viscous  and  source  term 
Jacobian  matrices  to  the  left  side  of  the  difference  scheme  is  required  to  produce  a  fully 
second-order  accurate  transient  solution. 

The  object  of  the  integration  scheme  is  to  produce  *.  This  is  actually  accomplished  at 
a  typical  point  as  described  in  what  follows.  Rewrite  Eq.  (24a)  as 

Offi,  =  RHS  (29) 

where  Bj,  B,  represent  the  factored  operators  on  the  left  side  of  Eq.  (24a)  and  RHS 
represents  everything  on  the  right  side  of  the  equal  sign.  Define  ij>*  such  that 

B^  =  0*  (30J 

Then 

B{0*  =  RHS  (31) 

and  the  algorithm  implementation  is  as  follows. 
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1 .  Given  0  at  all  points  at  level  n 

2.  Compute  RHS  and  elements  of  flj  in  Eq.  (31) 

3.  Solve  block  tridiagonal  system  to  get  <>* 

4.  Compute  elements  of  Q,  in  Eq.  (30) 

5.  Solve  Eq.  (30)  for  A<t> 

6.  Solve  for  =  0"  +  A<^ 

7.  Go  to  Step  2. 

3,4  CONSERVATIVE  PROPERTY 

It  was  reported  in  Ref.  21  that  exact  satisfaction  of  the  conservative  property  is  essential 
for  computation  of  flows  with  shock.  (It  is  unclear  to  this  author  that  such  a  requirement  is 
strictly  necessary  for  higher-order  differencing  schemes.)  This  requirement  means  that  a 
difference  scheme  must  numerically  satisfy  the  following  integral  condition  (divergence 
theorem) 


j  i  j  V  .  F  dV  =  §  n  •  F  dA  (32) 

V  SV 

where  F  is  the  flux  dyadic  represented  in  Cartesian  space  by  F  =  ei  +  ^.  For  a  2-D  problem 
in  jf  coordinates  the  discrete  analog  of  the  left  side  of  this  equation  is 

E  E  J  V  .  F  (33) 

j  k 


where  dV  =  Jd$  djj  and  df,  and  dvj  are  taken  to  be  unity  with  no  loss  of  generality.  The 
operator  JV  is  given  by 


JV  = 


Therefore, 
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JV  .  F 


de 


ae 


df 


df 

^dr, 


The  present  flux  difference  splitting  concept  modifies  this  to  the  following 


JV  .  F  =  Q^(y,ej  -  x^f^)  +  Q-(y,ej  -  x,f^) 


+  R-+(x/,  -  y^)  +  -  y^,) 

which  by  substitution  into  Eq.  (33)  gives  the  discrete  analog  of  the  left  side  of  Eq.  (32)  as 
f  s  [Q*(y,Sj  -  x,f,)  +  Q-(y,  Sj  - 
+  (xf  f,  -  yj  e,)  +  R-  (x^f,  - 

where  S  ^  indicates  a  double  summation  over  all  points  in  the  region  of  interest.  Now, 
inspection  of  Eq.  (32)  reveals  that  when  the  double  summation  in  Eq.  (34)  is  expanded  the 
result  must  be  that  all  terms  drop  out  except  terms  on  the  boundary  of  the  domain  of 
interest.  This  is  simply  a  fact  stated  by  the  divergence  theorem.  However,  in  practice,  this 
does  not  happen  exactly  unless  special  consideration  is  given  to  the  evaluation  of  the  metrics 
and  the  matrices  Q±,  R=.  This  special  consideration,  unfortunately,  depends  on  the 
differencing  formulae  selected  for  e^,  etc.  For  example,  consider  the  case  of  first-order 
spatial  differences  in  the  fluxes  e,  f.  Then  Eq.  (34)  becomes 

f  -  R-y^A^e) 

(35) 

-EE  (q"  W  +  Q-x,A/  -  RH-XjV/  -  R-x,A,f) 


where  terms  involving  e  and  f  have  been  separated.  An  expansion  of  the  double  summation 
shows  that  in  order  to  get  interior  point  cancellation,  the  metrics  x,,  yj,  y,  must  be 

centrally  differenced  and  the  products  Q+y,,  Q-y^.  etc.,  are  evaluated  with  two  point 
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averages  over  the  points  involved  with  the  associated  finite  difference.  For  example,  to 
compute  Q^y,  e,  the  product  Q^y,  is  evaluated  as 

(Q*!-,)  =  +  (Q^yOi-l]  0® 

This  assures  one  of  satisfying  global  conservation.  For  the  case  of  second-order  finite 

difference  representations  of  the  flux  derivatives  in  Eq.  (34),  no  method  has  been  found  of 
exactly  satisfying  the  global  conservation  property.  However  if  the  coefficients  Q'''y,.  Q"y,, 
etc.,  in  Eq.  (34)  are  evaluated  at  the  point  of  application  of  the  associated  second-order  flux 
difference  and  the  summation  indicated  is  carried  out,  then  it  can  be  shown  using  Taylor 
expansions  that  a  typical  interior  point  has  a  residual  which  reduces  to 

y  Af  +  higher  order  terms 

This  is  consistent  with  the  fact  that  the  second-order  difference  approximations  are  used 
at  each  point.  This  error  in  satisying  the  conservative  property  exactly  is  thought  to  be 
negligible  for  the  second-order  scheme,  although  further  tests  are  required  to  prove  or 
disprove  this  hypothesis. 

3.5  BOUNDARY  CONDITIONS 

This  section  describes  the  discretization  of  the  boundary  conditions  outlined  in  Section 
2.4.  The  general  philosophy  for  all  boundaries  is  to  incorporate  the  boundary  condition 
implicitly  and  then  impose  it  explicitly  after  the  integrated  solution  vector  is  determined. 
Refer  to  Fig.  5  for  the  following  discussion.  Consider  applying  the  factored  difference 
scheme  to  all  interior  points  in  the  computational  grid.  In  particular,  perform  the  step 
indicated  by  Eq.  (31)  along  each  tf  (or  k)  =  constant  line  excluding  k  =  1  and  k  =  K.  When 


Figure  5.  Typical  ri  =  constant  line. 
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the  operation  indicated  by  Eq.  (31)  is  expanded,  the  result  will  be  a  block  tridiagonal  system 
of  equations  to  solve  for  4>*  at  ail  points  from  j  =  2toj  =  J-  1.  This  system  has  the  form 

+  D2<I>2  +  U2(^3  =  RHSj 
+  U3^4  =  RHS3 


+  Dj-I0l-1  +  Uj-I0j  =  RHSj_i 

where  L,  D,  and  U  are  4  x  4  matrix  blocks  representing  lower,  diagonal,  and  upper  blocks, 
respectively.  This  system  involves  at  j  =  1  and  0*  at  j  =  J.  At  j  =  1,  or  boundary  D  of 
Fig.  2,  the  inflow  data  is  prescribed  and  never  changes,  implying  that  A0  =  0  here.  Since  <f>* 
is  proportional  to  A0  as  indicated  by  Eq.  (30),  then  0*  at  j  =  1  is  always  zero.  At  the  j  =  J 
boundary  (boundary  B  of  Fig.  2)  the  value  of  0j  is  needed.  In  the  present  work  0j  is  assumed 
to  depend  on  upstream  information  in  the  following  way 


0j  =  Of0j_,  +  00J. 


(37) 


If  (a,  ^  are  taken  to  be  (2,  -  1)  this  is  equivalent  to  enforcing  0jj  =  0  at  this  boundary,  (a, 
/S)  =  (1,0)  implies  that  0^  =  0  at  this  boundary.  Incorporation  of  these  conditions  at  j  =  1 
and  j  =  J  results  in  the  following  system  of  equations. 


Dj  Uj  0  0  0 

L3  Dj  Uj  0  0 


0  L4  D4  U, 


-J-3  ^J-3 


0  L 


J-2 


Uj-: 


0 


u, 


^J-2  *^3-2 

Lj_|  +  |6Uj_j  Dj-i  +  aUj. 


02 

4>*3 


RHS2 

RHS3 


0I-I  / 


rhSj_2 

RHSj.,/ 

(38) 


This  block  tridiagonal  system  is  solved  for  0^  through  0j_ ,  and  then  0j  is  computed  from 
Eq.  (37).  This  is  repeated  for  each  k  from  2  to  K  -  1. 
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The  approach  just  outlined  for  treatment  of  the  outflow  boundary  is  not  physically 
correct  even  from  a  signal  propagation  point  of  view.  This  is  because  subsonic  flow  will 
always  exist  in  the  boundary  layer  and,  technically  speaking,  some  mechanism  for  upstream 
influence  is  necessary.  However,  the  present  approach  does  work  for  the  class  of  problems 
considered  with  minimal  degradation  of  the  flow  upstream  of  the  outflow  boundary. 

Next,  a  solution  to  Eq.  (30)  for  A0  is  required  for  j  =  2  through  J.  Consider  a  typical  j 
=  constant  line  as  shown  in  Fig.  6.  Along  this  line,  expansion  of  Eq.  (30)  yields 


£.2^4^^  D2A^2  U2^03  —  02 

(39) 

LjA02  0^002  +  IJ2A04  =  03 

Lk-1^^K-2  +  Dk_iA0k_,  +  Uk^|A0k  =  0K-1  (40) 

This  system  involves  A0i,  and  A0k  which  are  wall  and  upper  boundary  increments  in  the 
solution  vector,  0.  The  treatment  of  A0|  is  addressed  first.  The  value  of  A0,  must  come 
from  the  wall  boundary  conditions  as  described  in  Section  2.4.  Since 

A01  =  (ap,  A(eu),  A(qv),  A(ee))^  (41) 

from  Eq.  (5a),  the  values  of  each  of  these  increments  at  the  wall  point  must  be  determined. 
The  no  slip  condition  requires  A(qu)  =  Afgv)  =  0.  The  density  increment,  Aq,  is  expressed 
in  terms  of  temperature  and  pressure  increments  from  the  equation  of  state  giving 

Aq - ^  Ap  -  -5-  AT  (42) 

^  RT  T 

The  normal  pressure  gradient  condition  is  used  to  obtain 

<4-)  = « = 


and 

A(.)  =  v(.)-n 

on 

where  fi  =  V7j/|V)j|  is  the  wall  normal  unit  vector.  So  since  V  =  Vf  +  Vij  fl/dij. 
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.  r,)  f  (Vv  .  V,,)  =  0 

iiii 

which  is  finite  differenced  to  yield  (see  Fig.  6) 

•  Vtj)j^4p  -  Apj_,J  +  y  (Vt;  •  V»;)j4Ap2  “  ^Ap  -  Ap,J  =  0  (43) 

Note  that  (his  equation  allows  expression  of  Ap  at  (he  wall  in  terms  of  Ap  at  points  on  the 
interior  adjacent  to  the  wall  and  the  first  upstream  wall  point. 


.1  -  1  .j 

Figure  h.  Typical  £  =  constant  line. 
The  temperature  condition  i.s  handled  as  follows.  Since 

o—  +  (I  -  <0(T  -  TJ  0 
^11 

time  differencing  yields 


a  — (AT)  +  (I  -  «)AT  =  0 
dn 

and 

>>AI  _  ^  «  Vt?  .1(41) 

f)n  l^’tl  fli  I  ^*7 1 


The  differenced  result  is 
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(4(AT)2  -  3(AT)  -  (AT)3)  -  AtJ 


<44) 


+  AT  =  0 


This  equation  allows  expression  of  AT  at  the  wall  in  terms  of  AT  at  points  on  the  interior 
adjacent  to  the  wall  and  the  first  upstream  wall  point.  Recall  that  a  =  0  for  a  fixed  wall 
temperature  and  and  a  =  1  for  an  adiabatic  wall  condition. 

The  quantity  A^|,  in  Eq.  (39)  must  be  related  to  the  wall  boundary  conditions  expressed 
by  Eqs.  (43)  and  (44)  along  with  the  no-slip  condition.  This  is  done  by  expressing  in 
terms  of  A0  at  interior  adjacent  points  and  possibly  A0  at  other  selected  points.  First 
recognize  that  the  fourth  element  of  <f>  is  total  energy  per  unit  volume  expressed  in  terms  of 
total  energy  per  unit  mass,  e,  as 


ge  = 


7  -  I 


^  [(eu)2  +  (ev)^] 


A  time  differencing  of  this  equation  yields 


A(ee)  =  — - (  “  ^  ^Ag  +  uA(gu)  +  vA(gv) 

7-1  V  2  / 


At  the  wall,  due  to  the  no-slip  condition, 

Ap  =  (7  -  I)A(ge) 


At  points  other  than  the 
Ap  =  (7  - 


wall,  Ap  is  expressed  in  terms  of  the  elements  of  A0  as 
O^Afge)  +  ^  ^  ^  ^Ag  -  uA(gu)  -  vA(gv)J 
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Let  the  vector  m  be  defined  as 


=  (7  -  1)[ 


+  v2  IT 

.  -  U,  -  V,  Ij 


Further  define  the  vector  J  as 


so  that 


Now  from  Eq.  (42)  it  is  clear  that 


Ap  =  m  •  A0 


e=  (i,o,o,o)T 


Ag  =  t  •  A<l> 


AT  =  Ap  —  —  Afl 

CR  Q 


AT  ~  -  m  •  A0  —  —  I  •  AA 

eR  e 


1  —  T  T 

n  =  - m  -  —  ( 

«R  Q 


Then 


AT  =  n  •  A0 


Substitution  of  Eqs.  (45b)  and  (46b)  into  Eqs.  (43)  and  (44)  gives 

am  .  A^  =  (V$  .  V»j)mj_,_,  •  A^j_,  ,  -  2iVj,  •  Vi/)mj2  •  Ai^j^j 


+  —  (Vtj  .  Vj,)mj_3  .  A0j_3 
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and 


(i>n  «  =  Q!(V{  •  1  ■  A<^j_i  1  -  2a(V’j*  VTf)nj_2  •  ^^j.2 


+  -7  (Vr>  •  Vi|)nj_3  .  A«^j_3 


(48) 


where 


a  =  (Vf  .  V,)  -  -|<V7,  .  Vt,) 


and 


CO  =  (1  —  a)  +  aa 

It  is  desirable  to  define  the  4  x  4  matrices  Z,  X,  W,  Y  such  that  the  boundary  condition  may 
be  expressed  as 


ZA^j  I  =  XA0j_|  I  +  WA0j  2  +  YA^j  3  (49) 

Equations  (47)  and  (48)  constitute  the  first  two  rows  of  Eq.  (49)  so  the  first  two  rows  of  the 
matrices  Z,  X,  W,  Y  may  be  written  down  by  inspection.  The  last  two  rows  of  Eq.  (49) 
simply  represent  the  no-slip  condition.  Let  the  four  elements  of  m  be  represented  by  m„  m2, 
m3,  and  m4  and  likewise  for  n.  Then  the  matrices  Z,  X,  W,  and  Y  are  given  by 


lami 

am2 

am3 

£UT1^^ 

conj 

con2 

<003 

0 

1 

0 

0 

\  0 

0 

1 

0 

(50a) 
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W  =  -  2(V,  .  V,)j  1 


Y  =(1/2)(V,  .  V,). , 


m2 

m3 

m4 

on2 

an3 

004 

0 

0 

0 

0 

0 

0 

f  mj 

m2 

m3 

an, 

002 

003 

0 

0 

0 

0 

0 

0 

m. 

m2 

m3 

m, 

an2 

003 

0 

0 

0 

0 

0 

0 

0  /  j,2  (50c) 


0  /  j,3  (50d) 


Now  define  X  -  Z  'X,  W*  =  Z“‘W,  Y’  =  Z  *Y.  Then  the  solution  increment  , 
becomes 

=  X*A0j_j  1  +  W*A0j2  +  Y*A<^j_3  (51) 

Thus,  with  this  result,  Eq.  (39)  is  modified  to  read 

(D2  +  L^W  )A^2  +  (U2  +  L2Y*)A03  =  —  L2X*A^j_]  I  (52) 
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Since  the  j  =  2  {  =  constant  line  is  computed  first,  the  value  of  ]  which  is  needed  on 
the  right-hand  side  of  Eq.  (52)  is  known  since  it  falls  on  the  inflow  boundary.  Each  $ 
=  constant  inversion  in  turn  requires  at  the  previously  computed  wall  point,  so  as 

long  as  the  sweep  proceeds  from  left  to  right  the  quantity  ,  in  Eq.  (52)  is  known. 

Next,  the  value  of  A^^jj  is  needed  in  Eq.  (40).  Two  approaches  are  considered  for  treating 
this  upper  boundary.  The  first  is  to  consider  the  upper  boundary  one  of  fixed  known  flow 
properties.  In  this  case,  A0jf  =  0  and  the  block  tridiagonal  system  of  equations  to  solve  at 
each  £  =  constant  line  from  j  =  2  to  j  =  J  is  given  by 


|D2  +  LjW*  U2  +  LjY*  0  0 

L3  D3  U3  0 


0 

L4  D4 

U4 

0 

0 

0 

Lk-3 

1 

Q 

Uk-3 

0 

0 

^K-2 

*^-2 

0 

0 

0 

^K-1 

Solution  of  this  system  yields  A^  through  A^^.j  for  the  given  j.  A^i^  is  zero  and  A«^3  is 
computed  from  Eq.  (St).  The  dependent  variable,  is  incremented  giving  the  new 
solution  on  the  given  j  =  constant  line.  Due  to  roundoff  errors,  this  solution  will  not  exactly 
satisfy  the  wall  boundary  conditions  even  though  these  conditions  are  implictly  satisfied.  For 
this  reason,  the  wall  conditions  are  explictly  enforced  with  the  equations 


Pj.i  = 


[ 


(V(  .  Vy)  pj_j  J  +  ^(Vt)  •  rjf)(pj_3  -  4  Pj  2) 


(V£  •  V,)  -  4  (Vti  .  Vtj) 


1 


(54) 


^  -r(Vf  .  -f  |(V,  .  T,)(Tj,3  -  4Tij)-|  ^ 

^j.l  “  “I  1  I  +  (1  “  Of)Tw 

L  (Vf  .  V,,)  -  4(V,.  Vi,)  J 


(55) 
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Uj_,  =  Vj ,  =  0  (56) 

Equations  (54)  -  (56)  provide  wall  information  to  be  used  to  compute  a  new  Even 
though  this  recomputed  wall  solution  is  only  slightly  different,  this  practice  of  reinforcing 
the  wall  boundary  conditions  explicitly  will  oftentimes  avoid  drifting  in  the  solution. 

The  second  approach  in  determining  is  to  actually  solve  the  governing  equations  on 
this  upper  boundary  since,  strictly  speaking,  the  first  approach  of  assuming  =  0  is 
physically  incorrect  from  a  signal  propagation  point  of  view.  The  first  thing  to  recognize  in 
this  context  is  that  must  be  determined  at  k  =  K,  which  requires  a  block  tridiagonal 
inversion  along  this  upper  boundary.  That  is,  Eq.  (38)  must  be  solved  along  k  =  K.  To 
present  a  clear  picture  of  the  approach  used,  it  is  necessary  to  begin  with  the  governing 
equations.  Define  the  vector  S  such  that 

C  _  i.  _  ^  I-  t. 

*  ”i  'y  di  dt)  "y  (57) 

The  governing  equation  is  then  written  [see  Eq.  (19)]  as 


+  S  =  0 

Premultiplication  of  this  equation  by  the  left  eigenvector  matrix  T  “ '  exposes  the  normalized 
eigenvalue  matrices  A*  on  the  ij  derivative  terms.  The  eigenvalues  associated  with  this 
direction  are 


>^1,2  =  ’ll  +  7j,U  -t-  V  =  +  Q  •  ^7 


Xj  =  Xi  +  c  ^ 
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X4  =  X,  -  c  VjJx  +  ’n] 

Since  Ht  =  -  g  •  Vri  where  g  is  the  speed  vector  of  the  point  in  question  (important  only  in 
moving  grid  problems),  X|  is  expressed  as 

5^1  =  (q  -  g)  •  =  Qr  • 

where  represents  the  fluid  velocity  relative  to  the  moving  boundary  point.  It  is  also  true 
that  Vi;  is  a'  vector  perpendicular  to  the  tj  =  constant  boundary  in  question.  In  this  case  Vi; 
points  out  the  computational  domain.  Therefore,  the  sign  of  X,,  indicates  whether  flow  is 
into  the  computational  domain  or  out  of  it  at  this  n  =  constant  boundary.  If  X,  >  0,  flow 
exits  the  computational  domain  through  the  upper  boundary  and  if  X,  <0  flow  enters  the 
domain  through  this  boundary.  Now  consider  X3. 

^  =  Qr  •  +  C  V»Jx  + 

or 

X3  =  V’lx  +  Vy  [Qr  •  »  +  c] 
or 

>^3  =  ^  [^RN  +  c]  =  C  +  7»y  (Mr^  +  1) 

r  Likewise, 

K  =  y/vl  +  [qrn  -  c]  =  cyjt}l  +  -  1) 

where  qjjj^  represents  the  relative  normal  velocity  through  the  upper  boundary.  Note  that 
qRjM  >  0  flow  exits  the  domain  and  qRj^  <  0  when  flow  enters  the  domain.  Thus  the 
sign  on  X3,  X4  is  determined  by  the  relative  normal  Mach  number  at  the  upper  boundary.  For 
the  problems  considered  in  the  present  work  the  relative  normal  Mach  number  is  a  low 
subsonic  value  so  that  X3  >  0  and  X4  <  0.  In  addition,  the  decision  was  made  to  treat  the 
upper  boundary  as  an  outflow  boundary  rather  than  an  inflow  boundary.  Thus  Xj  2  >  0 
which  gives  three  positive  eigenvalues  and  one  negative  eigenvalue.  The  compatibility 
equation  system  is 
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aci 


3{ 


(58) 


where  on  the  boundary  in  question 


a;  =  diag  (1.  1.  1.  0) 


and 


A,-  =  diag  (0,  0.  0,  1) 


The  ij-derivatives  associated  with  the  A,^ are  backward  difference  approximations  which  are 
obtainable  on  the  upper  boundary.  Therefore,  consistent  with  A^X  the  first  three  equations 
of  the  system  represented  by  Eq.  (58)  are  retained.  However,  the  fourth  equation  requires  a 
forward  ri  difference  on  the  A^terms  which  is  not  possible.  Therefore,  this  fourth  equation 
must  be  replaced  by  a  boundary  condition  on  this  upper  boundary.  The  condition  chosen  is 
to  enforce  a  “known”  static  pressure.  This  upper  boundary  pressure  is  not  really  precisely 
known  in  the  general  case  but  a  judicious  placement  of  the  upper  boundary  allows 
justification  for  the  use  of  this  condition.  The  result  is  that 

p  =  constant  in  time 

Thus,  Pt  =  0  and  this  condition  must  replace  the  fourth  compatibility  equation.  Since 


P  -  (7  -  1)  [♦.  -  J 


where  <t>i,  ^4  elements  of  the  solution  vector,  0,  p^  becomes 


[ 


P.  =  (y  -  1)1  <^>4.  -  —<f>2,  -  — <^31  + 


0| 


0) 


20f 


.] 
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or  P(  =  m  •  where  m  is  defined  by  Eq.  (45a).  Since  pj  =  +  i^p^  +  p,  the 

pressure  condition  is  expressed  as 


m  •  +  ft •  ^r,  =  ^  (5®) 

Equation  (59)  must  replace  the  fourth  compatibility  equation  of  the  system  represented  by 
Eq.  (58).  This  is  accomplished  by  first  defining  a  nonsingular  4x4  matrix  D“  S  which  has 
the  first  three  rows  identical  to  T“^  The  last  row  of  is  the  vector  m.  Also  let  T“* 
represent  T”‘  with  the  last  row  replaced  by  zeros.  With  these  definitions,  Eqs.  (58)  and  (59) 
are  expressed  as 

D-^r  +  +  D-‘Q-f,0j  +  D- V, 


T,-“  +  ^-w)  ^  ® 


+ 


/  de,  ^  dfi 

(’■-JT  ’'“a? 


)] 


=  0 


where  all  ij  derivatives  are  backward  in  sense.  Premultiplication  by  D  yields 

*r=  -  ® 

+  Q<£.-^  +  £,-|-)  *  ^  ] 

which  is  expressed  in  factored  finite  difference  form  as 

ft  Vjl  +  DT-**  Q+^£,VjA  + 


(60) 


+  |8AT^Q-f,A^  +  DT-*‘  Q-(fAA  +  j  j  x 
+  (36r^T/tV,I  +  Dr;i'  (ikV,A  +  T,yV^)J  ^ 
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-  +  +  DT;,i  Q  + 

+  Q-A^,Aj<^  +  DT,-*‘  Q-^Af^AfCi  +  AfyAjfj)  (61) 

+  Ai?,V^  +  DT,"#‘  (Ai?,V,ei  +  AjjyV^fj) J 

+  3A7^A{^5^»,  +  +  Aij^i^e,  +  A»fyi,fy^ 

where  the  same  approximations  used  in  arriving  at  Eq.  (24a)  for  interior  points  were  used 
here.  In  practice,  the  viscous  terms  in  Eq.  (61)  are  not  included  since  they  are  expected  to  be 
negligible  at  the  upper  boundary.  Equation  (61)  may  be  expressed  as 

Qfn,,A0  =  RHS 


or 

a^*  =  RHS 

This  latter  form  is  used  to  compute  <f>*  along  the  k  =  K,  i;  =  constant  boundary  as  indicated 
for  other  ij  -  constant  lines.  Then  since  the  operator  (1,  is  one-sided,  the  block  tridiagonal 
system  to  solve  for  A0  along  each  (  =  constant  line  is  written  by  appending  the  boundary 
equation  to  the  system  expressed  by  Eq.  (S3)  yielding 
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This  system  is  solved  for  through  A^^  £  =  constant  line  beginning  with  j  =  2 

and  ending  with  j  =  J.  The  wall  point  is  updated  as  discussed.  The  vector  is  modified  to 
assure  satisfaction  of  the  implicitly  enforced  pressure  condition.  This  is  done  by  simply 
accepting  the  first  three  elements  of  <^ic  and  recomputing  the  fourth  element  using  the  correct 
pressure  value.  If  this  step  is  not  performed  it  is  possible  for  the  upper  boundary  pressure  to 
deviate  from  its  desired  value  due  to  round  off  error.  This  is  particularly  true  when  a  large 
number  of  integration  steps  is  taken. 

3.6  INITIAL  CONDITIONS 

Before  the  governing  equations  are  integrated  using  the  techniques  described  in  the 
present  work,  the  computational  domain  must  be  initialized  to  some  starting  solution.  For 
the  class  of  problems  considered,  this  starting  solution  is  somewhat  arbitrary  since  only  the 
steady-state  solution  is  of  interest.  It  is  first  assumed  that  the  solution  on  the  inflow 
boundary  is  known  by  some  means.  This  could  be  data  supplied  by  another  computer  code, 
for  example.  Then  the  solution  on  each  $  =  constant  line,  in  turn,  is  set  equal  to  the  solution 
on  the  upper  boundary.  This  procedure  is  simple  and  reliable  for  the  class  of  problems 
considered. 


4.0  APPUCATIONS 

Two  computer  codes  were  written  using  the  procedures  outlined.  The  first  is  called 
VISCODE  and  the  second  is  VISCOD2.  VISCODE  imposes  a  fixed  boundary  condition  in 
all  dependent  variables  on  the  upper  boundary  of  the  computational  domain.  VISCODE2 
imposes  a  fixed  condition  only  on  static  pressure,  and  the  remaining  variables  are  computed 
on  the  upper  boundary.  VISCODE  was  applied  to  a  Couctte  flow  problem  and  V1SCOD2 
was  applied  to  a  flat  plate  laminar  boundary-layer  problem. 

4.1  COUETTE  FLOW 

Figure  7  depicts  a  schematic  of  a  Couette  flow  problem.  The  problem  tested  was 
characterized  by  the  following  parameters:  U„  =  100  meters  per  sec,  Po,  =  10*  Newtons  per 
square  meter,  T^  =  290.36  "K,  and  ymax  =  4  m.  A  coarse  5x5  grid  was  used  to  solve  this 
problem.  This  problem  has  an  exact  solution  when  coefficients  of  viscosity  and  thermal 
conductivity  are  constant.  The  exact  solutions  for  velocity  and  temperature  distributions  are 
given  by 


2k  “V 


ymax 


)(■  - 
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for  a  fixed  wall  temperature,  Tw.  and 

for  an  adiabatic  wall.  Figures  8  through  10  illustrate  a  comparison  of  the  numerical  solution 
versus  the  exact  solution  for  temperature  profiles  in  a  fixed  wall  temperature  environment. 
The  three  fixed  wall  temperatures  tested  were  T*  =  0,5  T^,  T*  =  T„,  and  T*  =  2T^. 
Figure  11  depicts  the  temperature  profile  comparison  for  the  adiabatic  wall  case.  Figure  12 
illustrates  the  velocity  profile  comparison. 


Figure  7.  Couette  flow  schematic. 


290.2  290.6  291.0  291.4 


Temperature,  T  (^K) 

Figure  8.  Temperature  profile,  Tw  =  T^ 
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This  test  problem  was  chosen  as  the  first  step  in  the  code  validation  process.  As  Figs.  8 
-  12  illustrate,  the  numerical  solution  is  right  on  target  for  this  problem. 


Temperature,  T  (^K) 

Figure  10.  Temperature  profile,  =  2  T„. 

4.2  BLASIUS  FLAT  PLATE  BOUNDARY  LAYER 

An  incompressible  laminar  boundary  layer  on  a  flat  plate  was  chosen  as  the  second  code 
validation  case.  See  Fig.  13  for  a  schematic.  This  problem  has  an  “exact”  solution 
formulated  in  terms  of  a  solution  to  a  third-order  ordinary  differential  equation.  This 
Blasius  solution  is  well  known  and  provides  an  additional  comparison  case  against  the 
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Temperature,  T  (OK) 

Figure  11.  Temperature  profile,  adiabatic  wall. 


Sym 

— ^  Exact 

O  Present  Calculation 


Figure  12.  Representative  velocity  profile. 
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Figure  13.  Flat  plate  problem  schematic. 

numerical  solution.  VISCOD2  was  used  to  compute  the  numerical  solutions.  Even  though 
the  exact  solution  used  is  for  an  incompressible  flow,  the  solution  to  the  compressible 
equations  should  compare  favorably  if  the  test  case  Mach  number  is  sufficiently  small.  The 
test  case  chosen  is  characterized  by  the  following  parameters:  Rcl  =  10^,  =  100  ft/ sec, 

=  500  ®R,  and  Pr  =  1  where  ReL  is  the  Reynolds  number  based  on  the  length  of  the  flat 
plate  and  Pr  is  the  Prandtl  number.  The  plate  wall  temperature  is  fixed  at  600’R.  These 
conditions  translate  to  a  Mach  number  of  0.09,  which  is  sufficiently  small  to  justify  the 
comparison.  The  exact  solution  for  the  temperature  profile  as  a  function  of  local  velocity  for 
the  problem  described  is  given  by 


T(u)  =  (l  -  +  (-JL-)t.  +  ^1  - 


where  Cp  is  the  specific  heat  at  constant  pressure  of  air. 

A  comparison  of  known  versus  numerical  solutions  is  found  in  Figs.  14  and  15.  Figure  14 
depicts  the  velocity  profile  comparison  at  a  downstream  station.  Two  numerical  solutions 
are  presented  in  this  figure.  One  is  for  the  case  of  a  uniformly  spaced  grid  in  the  y-direction 
and  the  other  is  for  the  case  of  a  grid  clustered  toward  the  plate  to  the  extent  that  the  first 
point  off  the  wall  is  at  y  =  0.0005  ft.  In  both  cases,  30  points  in  the  y  direction  were  used. 
Figure  15  illustrates  the  comparison  for  the  boundary-layer  temperature  profile  at  a 
downstream  station.  All  comparisons  show  good  agreement.  As  a  result,  this  part  of  the 
code  validation  is  considered  to  be  successful. 
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Temperature,  T  (°K) 

Figure  15.  Temperature  profile  at  exit  plane. 
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5.0  CONCLUSIONS  AND  FUTURE  WORK 

A  two-dimensional  or  axisymmetric  Navier-Stokes  code  was  developed.  The  code  is 
formulated  using  a  factored  delta  form.  All  in  viscid  terms  are  treated  using  a  flux  difference 
splitting  concept.  The  viscous  terms  are  centrally  differenced.  A  fully  general  geometric 
description  is  employed  and  all  metric  terms  are  retained,  including  those  pertaining  to  a 
moving  grid.  The  algorithm  as  presently  constructed  is  first  order  in  time  and  either  first  or 
second  order  in  space.  It  can  be  made  seond  order  in  time  by  the  inclusion  of  the  viscous  and 
source  term  Jacobians  in  the  operators  on  the  factored  side  of  the  scheme.  The  code  shows 
very  promising  stability  properties  and  convergence  rate.  Two  codes  were  actually 
developed.  One  imposes  a  fixed  strategy  for  all  variables  on  the  upper  boundary  of  the 
computational  domain.  The  other  imposes  a  fixed  strategy  only  on  the  static  pressure  at  this 
upper  boundary.  Validation  runs  were  made  on  a  Couette  flow  problem  and  a  Blasius  flat 
plate  boundary  layer.  Results  were  good.  The  code  must  yet  be  validated  on  axisymmetric 
problems  and  problems  with  shocks  such  as  a  simple  ramp  induced  shock-boundary  layer 
interaction.  A  turbulence  model  is  currently  being  incorporated  to  extend  the  current 
capability  to  turbulent  flows.  Future  work  should  also  include  a  conclusive  study  of  the 
importance  of  the  conservative  property  and  possibly  development  of  an  adaptive  grid 
scheme  for  the  class  of  problems  of  interest. 
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NOMENCLATURE 

A  Jacobian  matrix  associated  with  x-direction  inviscid  flux 

a  Coefflcient  used  in  application  of  wall  boundary  condition 

A*  Jacobian  matrix  associated  with  ^-direction  inviscid  flux 

B  Jacobian  matrix  associated  with  y-direction  inviscid  flux 

B*  Jacobian  matrix  associated  with  ij-direction  inviscid  flux 

C  Jacobian  matrix  associated  with  x-direction  viscous  flux 
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c  Sonic  velocity 

Also  a  solution  dependent  coefficient  used  to  explain  viscous  term  differencing 

Cp  Speciflc  heat  at  constant  pressure 

D  Jacobian  matrix  associated  with  y-direction  viscous  flux 

Also  diagonal  matrix  in  block  tri- diagonal  system 
Also  inverse  of  boundary  condition  altered  matrix  of  left  eigenvectors 

D  ■  *  Boundary  condition  altered  matrix  of  left  eigenvectors 

e  Specific  total  energy 

Also  y-direction  flux  vector  in  cylindrical  frame 

e  y-direction  flux  vector  in  Cartesian  frame 

Ci  Inviscid  part  of  y-direction  flux  vector  in  cylindrical  frame 

Ci  Inviscid  part  of  y-direction  flux  vector  in  Cartesian  frame 

Viscous  part  of  y-direction  flux  vector  in  cylindrical  frame 

Cv  Viscous  part  of  y-direction  flux  vector  in  Cartesian  frame 

F  Dyadic  flux  function  for  Navier-Stokes  equations 

f  x-direction  flux  vector  in  cylindrical  frame 

f  x-direction  flux  vector  in  Cartesian  frame 

f]  Inviscid  part  of  x-direction  flux  vector  in  cylindrical  frame 

fi  Inviscid  part  of  x-direction  flux  vector  in  Cartesian  frame 

fv  Viscous  part  of  x-direction  flux  vector  in  cylindrical  frame 

fv  Viscous  part  of  x-direction  flux  vector  in  Cartesian  frame 

G  Grid  clustering  function 
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g  Geometry  related  coefficient  used  to  explain  viscous  term  differencing 

g  Grid  point  speed  vector 

h  Source  term  in  cylindrical  frame 

Hj  Jacobian  matrix  associated  with  inviscid  source 

hi  Inviscid  part  of  source  term  in  cyindrical  frame 

Hy  Jacobian  matrix  associated  with  viscous  source 

hv  Viscous  part  of  source  term  in  cylindrical  frame 

I  Identity  matrix 

i  x-direction  unit  vector 

J  Jacobian  of  coordinate  transformation 

Also  number  of  grid  points  in  the  {-direction 

j  {-direction  index 

A 

j  y-direction  unit  vector 

K  Number  of  grid  points  in  the  ij-direction 

k  Coefficient  of  thermal  conductivity 

Also  i]-direction  index 

L  Lower  off-diagonal  matrix  in  block  tri-diagonal  system 

f  Vector  used  in  application  of  body  boundary  condition 

m  Vector  used  in  application  of  body  boundary  condition 

^R|si  Relative  normal  Mach  number  used  at  upper  boundary 

n  Direction  normal  to  a  surface 

Unit  vector  in  direction  of  n 
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n  Vector  used  in  application  of  body  boundary  condition 

P  {-forcing  function  in  grid  generation  scheme 

p  Static  pressure 

Pr  Prandti  Number 

Q  T^forcing  function  in  grid  generation  scheme 

q  Cylindrical  velocity  vector 

q  Cartesian  velocity  vector 

qR  Relative  velocity  vector 

Qrk  Relative  normal  velocity 

Q*  Rotation  matrices  used  to  control  differencing  of  {-direction  in  viscid  terms 

R  Gas  constant 

r  Position  vector 

R^  Rotation  matrices  used  to  control  differencing  of  i}-direction  inviscid  terms 

RHS  Right  hand  side  of  integration  algorithm 

S  Solution  dependent  quantity  used  to  explain  viscous  term  differencing;  total  arc 

length  along  f  =  constant  curve;  also  represents  grouping  of  source  and 
viscous  terms 

s  Arc  length  along  grid  lines 

T  Temperature 

t  Time 

Matrix  of  right  eigenvectors  associated  with  {-direction 
^  Matrix  of  left  eigenvectors  associated  with  {-direction 
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T,  Matrix  of  right  eigenvectors  associated  with  T^direction 

^  Matrix  of  left  eigenvectors  associated  with  i;-direction 

U  Upper  off-diagonal  matrix  in  block  tri-diagonal  system 

u  x-direction  cylindrical  velocity 

u  x-direction  Cartesian  velocity 

Free-stream  velocity 

V  y-direction  cylindrical  velocity 

V  y-direction  Cartesian  velocity 

W  Matrix  used  in  application  of  body  boundary  conditions 

W*  Matrix  used  in  application  of  body  boundary  conditions 

X  Matrix  used  in  application  of  body  boundary  conditions 

X  Independent  Cartesian  coordinate 

X*  Matrbc  used  in  application  of  body  boundary  conditions 

V  Matrix  used  in  application  of  body  boundary  conditions 

y  Independent  Cartesian  coordinate 

Y*  Matrix  used  in  application  of  body  boundary  conditions 

Z  Matrix  used  in  application  of  body  boundary  conditions 

z  Independent  Cartesian  coordinate 

a  Coefficient  used  in  grid  equations 

Also  used  in  extrapolated  outflow  boundary  condition 
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a  Used  as  switch  to  choose  temperature  boundary  condition 

13  Coefficient  used  in  grid  generation 

Also  a  parameter  used  in  clustering  the  grid 

Also  used  in  extrapolated  outflow  boundary  condition 

Also  used  to  switch  order  of  time  differencing 

^  Related  to  0  as  used  to  switch  order  of  time  differencing 

r  +  General  backward  difference  operator 

r-  General  forward  difference  operator 

7  Coefficient  used  in  grid  generation 

Also  ratio  of  specific  heats 

A,  V  Forward  and  backward  difference  respectively 

d  Half  point  central  difference  operator 

f  Transformed  independent  variable 

71  Transformed  independent  variable 

^  Normalized  transformed  independent  variable 

A  Diagonal  eigenvalue  matrix 

A+  Diagonal  eigenvalue  matrix  of  positive  eigenvalues 

A"  Diagonal  eigenvalue  matrix  of  negative  eigenvalues 

A*  Normalized  diagonal  eigenvalue  matrix  with  positive  or  negative  elements 

X  Second  coefficient  of  viscosity 

Also  eigenvalue  of  Jacobian  matrix 

/i  First  coefficient  of  viscosity 

{  Transformed  independent  variable 
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Q  Fluid  density 

£  Indicates  discrete  summation 

Also  used  in  extrapolated  outflow  boundary  condition 

T  Temporal  coordinate 

A  part  of  forcing  function  used  in  grid  generation 

^  Cylindrical  dependent  variable  vector 

0  Cartesian  dependent  variable  vector 

Intermediate  dependent  variable  vector 

\ff  A  part  of  forcing  function  used  in  grid  generation 

17-direction  factored  implicit  operator 

17-direction  factored  implicit  operator 

w  Coefficient  used  in  application  of  wall  boundary  condition 

Subscripts 

i  Inviscid 

j,k  Grid  point  indices 

t,x,y,7,{,i7  Indicates  partial  differentiation  in  associated  direction 

Also  indicates  simple  association  with  a  particular  direction 

v  Viscous 

w  Wall 
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Superscripts 

n  Temporal  index 

T  Transpose 

+  Associated  with  positive  eigenvalues 

Associated  with  negative  eigenvalues 

Other 

V  Gradient  operator 

Also  backward  finite  difference 

Laplacian  operator 
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